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Abstract 

We solve the equations of motion of a one-dimensional planar Heisenberg (or Vaks- 
Larkin) model consisting of a system of interacting macro-spins aligned along a ring. 
Each spin has imit length and is described by its angle with respect to the rotational 
axis. The orientation of the spins can vary in time due to random forcing and spin- 
spin interaction. We statistically describe the behaviour of the sum of all spins for 
different parameters. The term "domino model" in the title refers to the interaction 
among the spins. 

We compare the model results with geomagnetic field reversals and find strikingly 

similar behaviour. The aggregate of all spins keeps the same direction for a long time 
and, once in a while, begins flipping to change the orientation by almost 180 degrees 
(mimicking a geomagnetic reversal) or to move back to the original direction (mimick- 
ing an excursion). Most of the time the spins are aligned or anti-aligned and deviate 
only slightly with respect to the rotational axis (mimicking the secular variation of the 
geomagnetic pole with respect to the geographic pole). Reversals are fast compared 
to the times in between and they occur at random times, both in the model and in 
the case of the Earth's magnetic field. 

Keywords: Ising-Heisenberg model - geomagnetic field - reversals - statistical me- 
chanics 

1 Introduction 

One of the most remarkable phenomena of geomagnetism is that the Earth has reversed the 
polarity of its almost dipolar magnetic field many times in the past at irregular intervals 
(e.g., Jacobs, 1994; Merrill et al., 1996). Similar reversals have also been observed in 
turbulent dynamo experiments (Berhanu et al., 2007) and in simulations of the geodynamo 
(Glatzmaier and Roberts, 1995), but the cause of the reversals has yet eluded a convincing 
explanation (Amit et al., 2010). 
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The magnetic field of the Earth originates from dynamo action in the hquid outer 
core (e.g., Roberts, 2007). Numerical simulations of the geodynamo successfully repro- 
duce many features of the magnetic field of the Earth including stochastic reversals (e.g., 
Christensen and Wicht, 2007). Depending on the importance of the inertial forces rel- 
ative to the rotational forces, dynamos with either a dominant axial dipole or with a 
small-scale multipolar magnetic field are found (Kutzner and Christensen, 2002). The 
transition from dipolar to multipolar dynamos takes place at a local Rossby number of 
approximately 0.12 (Christensen and Aubert, 2006). The Earth lies close to the transition 
between both types (Olson and Christensen, 2006), which may explain why the dipole 
undergoes sporadic reversals. 

Convection in fast rotating bodies like planetary liquid interiors organises itself in the 
form of convective columns that encircle the inner core tangent cylinder and are aligned 
in the z-direction parallel to the rotation axis. The flow becomes quasi two-dimensional 
(geostrophic), minimising any variation in z-direction (Proudman- Taylor theorem). Cy- 
clonic and anticyclonic columns, rotating in the same and in the opposite direction as the 
planet, respectively, alternate in azimuthal direction (Busse, 1975). In strongly driven 
dynamos the number of columns increases. Secondary flows are directed away from the 
equatorial plane in anti-cyclonic columns, but towards the equatorial plane in cyclonic 
columns. The helicity, a key ingredient in the dynamo process, has therefore predomi- 
nantly one sign in the northern hemisphere and the opposite sign in the southern hemi- 
sphere. The so-called a^-dynamo mechanism, described by Kageyama and Sato (1997) 
and Olson et al. (1999), can be thought of as a process where the helicity associated to 
each convective column produces its own magnetic field. The alignment with the rotation 
axis and the organised helicity guarantees that the sum of these individual contributions 
add up to form the dominant axial dipole field. 

In the numerical simulations. Earth-like reversals are rare and only happen once the 
driving is large enough, whereby the fields arc very complex. There have therefore only 
been a few attempts to unravel the reason for the polarity switches (Glatzmaicr and 
Roberts, 1995; Kida and Kitauchi, 1998; Glatzmaier ct al., 1999; Sarson and Jones, 1999; 
Kutzner and Christensen, 2004; Wicht and Olson, 2004; Wicht, 2005; Takahashi et al., 
2005; Bouligand et al., 2005; Aubert et al., 2008; Wicht et al., 2009; Amit et al., 2010). 
Typically these simulations show that reversals go along with a breaking of the north- 
south (equatorial) symmetry in the fiow of the aligned fluid columns. Aubert et al. (2008) 
identified flow upwellings as the cause of held reversals. These features rise from the 
inner-core boundary and produce inverse magnetic fleld. To some extend, the upwellings 
behave like tilted convective columns, at least what their role in the dynamo mechanism 
is concerned. 

The reversal sequences in paleomagnetic data and in dynamo models have been anal- 
ysed for their statistical properties. First estimates that geomagnetic reversals obey a 
Poissonian process where all reversals are independent of each other have been shown to 
fall too short. The statistical reversal rate has likely been changing over time due to the 
varying heat flux through the core- mantle boundary (Glatzmaier et al., 1999; Constable, 
2000; Driscoll and Olson, 2009). The reversal sequence also suggests that the process 
may have a short and a long term memory, leading to changes in the statistical behaviour 
and the characterising distribution function of the times between reversals (Jonkers, 2003; 
Ryan and Sarson, 2007). Similar analysis for fully 3D numerical dynamo simulations are 
rare because it is very costly to compile a large number of reversals. The analysis by Wicht 
et al. (2009) and Driscoll and Olson (2009) indicate that the numerical simulations may 
follow a similar reversal statistics as the paleomagnetic record. 
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Simple parameterised models allow for a large number of reversals so that a statistical 
analysis becomes more meaningful. A famous example is the two-disk dynamo of Rikitake 
(1958), which exhibits sporadic reversals, but also a cyclic variation of the dipole moment 
during stable polarity periods. The extensions to N coupled disks by Shimizu and Honkura 
(1985) and Ito (1988) improved on the latter weak point. Hoyng ct al. (2001) considered a 
mean-field dynamo model with stochastic fluctuations of the induction effect; these lead to 
oscillations of the dipole field amplitude in a bistable potential with minima representing 
normal and reversed polarity and occasional jumps between them (Schmitt et al., 2001). 

Here we study another class of simplified models, an Ising-Heisenberg model of in- 
teracting magnetic spins. Ising-like models have been used in molecular dynamics and 
statistical mechanics for describing, for instance, phase transitions in ferromagnetism, for 
modelling spin glasses and for pattern recognition in neural networks (e.g., Greiner et al., 
1995). Coupled spin models of Ising type, where the individual spins can assume two 
scalar states +1 or —1 and interact with each other after certain rules, have also been 
suggested for describing geomagnetic polarity reversals and their statistics (Mazaud and 
Laj, 1989; Seki and Ito, 1993; Dias ct al., 2008). We analyse a planar Heisenberg model 
consisting of a system of vectorial spins aligned along a ring. Each spin has unit length and 
is described by its angle with respect to the rotational axis, i.e., each spin has one degree 
of freedom. The orientation of the spins can vary in time due to spin-spin interaction and 
random forcing. The consecutive interaction of adjacent spins is described as "domino 
model" . We consider the time dependence of the average orientation of all spins, which 
exhibits a similar behaviour as the geomagnetic reversal record. 

This sort of models are often classified according to spatial dimensionality and number 
of components of spin vectors. Our model is one-dimensional (i.e., spatial dimensionality 
of the lattice is one) and vectors (spins) are two-dimensional (i.e., they are contained 
within a plane). Thus, the domino model is a one-dimensional XY model, also referred 
to as the plane rotator model or the Vaks-Larkin model (Vaks and Larkin, 1966). Spin 
vectors in Ising and Heisenberg models are ID and 3D, respectively. A clear classification 
can be found in the classical textbook by Stanley (1987). 

The spins in our model might be associated with the convection columns, whose elec- 
tromagnetic induction generates elementary dipoles. The tendency of the spins to be 
aligned with the rotation axis is a consequence of the Proudman- Taylor theorem, and the 
time variation of the spins is a measure of the vigor of convection and of the sporadic up- 
wellings. The domino model may be thought of to mirror the full dynamo in a simplistic 
manner. One should, however, be careful not to take the analogies too far. The domino 
model is, after all, a toy model for a rather complex process. 

The structure of the paper is as follows. In Sect. 2 the model is described. The results 
and the statistical analysis of our model are presented in Sect. 3. The influence of the 
various parameters as well as alternative model descriptions are given in Sect. 4. In Sect. 5 
the results are compared with geomagnetic data, and in Sect. 6 we give our conclusions. 

2 The domino model 

We consider a system of N macro-spins aligned along a ring and interacting pairwise like 
in a one-dimensional Vaks-Larkin model. The spins are embedded in a uniformly rotating 
medium and we take ft = (0, 1) as the unit vector along the rotational axis. Each spin 
Si, i = 1,...,N has unit length and is described by its angle 9i with respect to the 
rotational axis, such that Si = {sin 6i, cos 9i). The orientation of the spins can vary in 
time due to random forcing and spin-spin interaction (Fig. 1). 
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Figure 1: Sketch of the domino model. 



The kinetic and the potential energy K{t) and P{t) of the system are 
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where i + 1 = 1 when i = N. Here 7 is a parameter characterising the tendency of 
the spins to be aligned with the rotation axis, while A is a parameter characterising the 
spin-spin interaction. The scalar product to the square in the 7-term ensures that there 
is no preferred polarity. The interaction is such that each spin interacts with the two 
neighbouring spins: spin 2 interacts with spins 1 and 3, spin 3 with spins 2 and 4 and so 
on. Spin interacts with spins — 1 and 1, i.e. we are considering periodic boundary 
conditions. 

The Lagragian for the system is C = K — P. We set up a Langevin-type equation as 
follows 

d f dC\ dC A / .N . eXi 



Kei{t) + 



(3) 



where the term —K9i{t) describes friction and the term exi/ y/r is a random force acting 
on each spin. The parameters k and e characterise the strengths of the friction and the 
random forcing, respectively. Finally, Xi is a Gaussian-distributed random number with 
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zero mean and unit variance associated to each spin, which is updated each correlation 
time T. 

Inserting the expressions for the kinetic and potential energy, (1) and (2), into equation 
(3) yields 

Oi — 27 cos Oi sin 9i + \ [cos 9i (sin + sin ^i+i) — sin Oi (cos di-\ + cos ^i+i)] 

+ /t^i-^ = 0, i = l,...,Ar (4) 

with ^0 = and ^at+i = 9i. 

We integrate the equations of motion (4) forward in time with a 4th-order Runge-Kutta 
scheme starting from a random orientation of the spins between and 27r. A standard set 
of parameters N , 7, A, k, e and r is considered in Sect. 3, while the parameter dependence 
of the mean time between reversals is studied in Sect. 4.1. In Sect. 4.2 we also slightly 
alter the model using different alternatives of the 7-term, the A-term and the forcing term. 

As main output we consider the cumulative orientation of all spins and define 

N N 

i=l 1=1 
as the resulting total axial magnetic moment or "magnetisation" . 



3 Results and statistical analysis of a typical model 

The parameters of our standard model are AT = 8, 7 = —1, A = —2, k = 0.1, e = 0.4 
and T = 0.01. The integration time step was At = 0.01, the total number of time steps 
was 30 000 000, and every 10th time step is outputted. The run comprises a total of 824 
reversals, i.e. the mean time between reversals is 364. Identifying this time with the mean 
time between reversals in the case of the Earth, which is 300 kyr, the whole nm with a 
time of 300 000 spans approximately 250 Myr. In Fig. 2 the first tenth of the full run is 
displayed. 
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Figure 2: The magnetisation of the standard run as a function of time. At the top the 
times of all zero-crossings are indicated (all reversals), while at the bottom only those 
where a central band of M = [—0.5,0.5] is crossed are displayed (true reversals). 

The statistical analysis is done using the whole run. The power spectrum is shown in 
Fig. 3a. Over a large range comprising most of the reversals, the spectrum follows a power 
law with an exponent of about —1.7. The spectrum for small frequencies or long polarity 
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Figure 3: Left (a): Log-log plot of the power spectrum of the magnetisation. Right (b): 
Normalised distribution of the magnetisation. 
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Figure 4: Left (a): Log- log plot of the number of reversals as a function of the times 
between reversals. Right (b): Normalised distribution of the times between reversals for 
short duration chrons. 



chrons (i.e., epochs of one polarity) is flatter, while the steeper decrease at high frequencies 
comes from the fast variations between reversals. The distribution of the magnetisation 
peaks near ±1 with a wide and deep valley between them (Fig. 3b). This reflects that the 
flipping or reversal times arc short events compared to the average diiration time between 
them and that the spins are most of the time closely aligned with the rotation axis. 

The distribution of the duration of long chrons follows a power law with an exponent 
of approximately —1.5 (Fig. 4a), while the short chrons are approximately exponentially 
distributed (Fig. 4b). 

In Fig. 5 the details of one single reversal at t ~ 3450 is shown. The reversal is triggered 
by a large fluctuation of spin 4 which is successively transferred to the neighbouring spins, 
which fluctuate until finally all spins reverse their polarity. We call our model "domino 
model" because of this consecutive interaction of neighbouring spins. The duration of 
the full reversal depends on how fast the original fluctuation is transferred to all other 
spins. The reversal here lasts about 10 time units, which corresponds to 8000 yr, roughly 
the time it takes for the geomagnetic field to perform a reversal. Shortly afterwards, at 
t ~ 3470, some spins again show large variations and even reverse for a short time, but 
they fail to transfer this to all the spins. The total magnetisation shows an excursion. 

Examples of a true reversal, an aborted reversal and two excursions are shown in Fig. 6. 
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Figure 5: The total magnetisation of all spins, cos6i/N, (top panel) and the mag- 

netisation of the individual spins, cos 9i, before and after a reversal at i 3450. 
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Figure 6: Examples of a true reversal (top panel), an aborted reversal (middle panel) and 
two excursions (bottom panel). 



4 Parameter study 

The model contains a number of free parameters. The main quantity which depends 
sensitively on their values is the frequency of reversals. The results, described in Sect. 4.1, 
are based on runs with at least three hundred reversals. The slope of the power spectrum 
depends only weakly on the model parameters. The distribution of the magnetisation 
is also similar to the model of Sect. 3. When there are very many reversals, the valley 
between the two stable states is less wide and deep. 
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Figure 7: The mean time between reversals depending on the number of spins N. The 
symbol for the standard model of Sect. 3 is plotted larger. The values of all other param- 
eters are as in our standard model. 
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Figure 8: The mean time between reversals depending on (a) 7, (b) A, (c) k, and (d) e. 
The symbol for the standard model of Sect. 3 is plotted larger. All other parameters are 
as in our standard model. 



4.1 Mean time between reversals 

The dependence of the mean time between reversals on the mimbcr of spins is shown in 
Fig. 7. In the case of a few spins this time steeply increases with the number, but it seems 
to saturate for a large number of spins. 

The parameter 7 measures the tendency of the spins to be aligned with the rotation 
axis. Large negative values stabilise the orientation and lead to fewer reversals (Fig. 8a). 



8 



A value of 7 = still behaves like a bistable oscillator, but with many reversals and an 
almost flooded valley between the stable states in the distribution function. Positive values 
of 7 do not lead to a stable magnetisation, but to oscillations around M = 0. 

The influence of the spin-spin interaction parameter A is similar. Large negative values 
stabilise (Fig. 8b). For A > the system randomly oscillates around M = 0. 

Increased friction, described by larger values of k, quite naturally stabilises (Fig. 8c), 
while increased random forcing, described by larger values of e, destabilises the system, 
leading to shorter chrons and more frequent reversals (Fig. 8d). 

4.2 Alternative model descriptions 

To test the robustness of the domino model we also slightly altered its detailed description. 

Instead of the 7-term in Eq. (2) proportional to J2{^ ' •^i)'^ have considered a 
term proportional to ^ |fi • Si\. Since all \ fl • Si\ < 1 and thus ^ |fJ • Si\ > X](fi • Si)'^ 
very similar results were obtained for just somewhat smaller (absolute) values of 7. For a 
7-term proportional to ^(il ■ Sj) no reversals were observed. 

Furthermore, instead of additive forcing with Gaussian noise, described by the last term 
in Eq. (3), we also investigated white noise. This results in no much change, except for 
slightly less frequent reversals. We also studied multiplicative forcing proportional to cos di; 
this results in much fewer reversals. Multiplicative forcing proportional to mod(6'i,27r) 
leads to random oscillations about 0. 

As an alternative to the local interaction with neighbouring spins only, described by 
the A-term in Eq. (2), we also considered a global or "mean-field" interaction with all 
other spins, described by {2X/N) ^^^^('S'i • Sj). A normalisation factor of 2/N is included 
in order to compare with the standard interaction with just the adjacent neighbours. 
The mean-field model results in less frequent reversals, but shows otherwise qualitatively 
similar behaviour. 

It is interesting to observe that in globally coupled models a qualitatively similar 
behaviour to the domino model described above is also found even without noise and 
without friction. The resulting system is conservative and only two parameters, 7 and A, 
are left. We shall study this system in detail in a forthcoming paper (Nakamichi et al., 
2011). 

5 Comparison with geomagnetic data 

We use the geomagnetic polarity time scale of Cande and Kent (1992, 1995) and Ogg 
(1995), which covers the past 166 Myr and comprises 332 reversals. Assigning a magneti- 
sation of -|-1 for chrons of normal polarity and —1 for crons of reversed polarity we derive 
Fig. 9, which is an analogue to Fig. 2. Reversals occurred at irregular intervals of 10^ to 
10'' yr. The mean time between reversals is approximately 300 kyr, whereas reversals are 
fast events lasting a few kyr. The reversal frequency has considerably decreased towards 
and increased away from the Cretaceous superchron which lasted from 118 to 83 Myr BP 
(Constable, 2000). 

The non-stationarity of the geomagnetic reversal record is not present in the domino 
model of Sect. 3, but could be easily accounted for by a gradual change of the model 
parameters with time (Sect. 4). It presents some difficulties in the statistical analysis of 
the record. The cumulative distribution of polarity chrons follows a power law with an 
exponent of —1.5 (Fig. 10a). Polarity intervals of a duration shorter than 1 Myr, which 
make up the vast majority of all intervals, follow an exponential or Poissonian distribution 
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Figure 9: Geomagnetic reversal record from present to 166 Mio yr BP. The record com- 
prises 332 reversals. A magnetisation of +1 is assigned for chrons with normal polarity 
and of —1 for chrons with reversed polarity. 



with a mean of 300 kyr (Fig. 10b). Ryan and Sarson (2007) find that the full set of polarity 
intervals is better fitted by lognormal and loglogistic distributions rather than Poisson and 
gamma distributions (Constable, 2000). The power spectrum of the geomagnetic record 
follows power laws with an exponent of about —0.6 for chrons longer than about 3 Myr 
and an exponent of about —1.9 for chrons of shorter duration (Fig. 10c). 

As a measure of the short-term variability of the geomagnetic field at times between 
reversals we analyse the virtual axial dipole moment (VADM) of the SINT-2000 data 
set (Fig. 11a) (Valet et al., 2005). The power spectrum of these fluctuations with a 
characteristic power index of —3 is displayed in Fig. lib, while the distribution of the 
VADM is given in Fig. 11c. For a comparison of the distribution derived from the shorter 
SINT-800 data set of Guyodo and Valet (1999), see Hoyng et al. (2002). 

When comparing the geomagnetic data with the behaviour of our reference model il- 
lustrated in Figs. 3 and 4 striking similarities become apparent. The fits to the chron 
durations suggest similar power law exponents of approximately —1.5 for the paleomag- 
netic sequence and for our simple model. In the latter we have disregarded the short 
chrons, which seem to represent brief statistical ventures into the other polarity because 
they follow a different behaviour. These may be identified with paleomagnetic excursions. 

The power spectrum of the geomagnetic reversal record (Fig. 10c) only refers to polarity 
epochs and therefore does not contain the high frequency contributions in our model 
(Fig. 3c). The low- frequency range can be interpreted as the background variation in the 
reversal frequency. This leaves us with comparing the mid-frequency spectrum with a 
slope of —1.9 for the paleomagnetic data and —1.7 for our model. The high frequency 
part of the model can be compared with the SINT data analysis, which yield a slope of —3 
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Figure 10: Top left (a): Log- log plot of the number of reversals as a function of times 
between reversals of the geomagnetic reversal record shown in Fig. 9. Bottom left (b): 
Normalised distribution of the times between reversals for short duration chrons of less 
than 1 Myr. These make up 298 out of the total of 332 reversals. The bin size here is 
four times the interval size. The dashed line is the expected probability density function 
in the case of a Poissonian process with a mean polarity residence time of 300 kyr. Right 
(c): Power spectrum of the geomagnetic reversal record. 

compared to —6 in the model. Not surprisingly, this discrepancy suggests that our model 
does a good job in replicating the reversal dynamics, but not in the details of the secular 
variation. A detailed discussion of the power spectrum of reversals as well as intensities is 
presented in Constable and Johnson (2005). 

6 Discussion and conclusions 

Our simple domino model of interacting magnetic spins reproduces the qualitative fea- 
tures of geomagnetic polarity reversals remarkably well. The orientation of the aggregate 
of all spins is most of the time nearly aligned or anti-aligned and deviates only slightly 
with respect to the rotational axis. Once in a while, at sporadic times, it starts flip- 
ping to ultimately change the orientation by almost 180 degrees or to move back to the 
original direction. The model thus mimics sporadic reversals of polarity, excursions and 
secular variation of the geomagnetic field. The power spectrum derived from the paleo- 
magnetic reversal records as well as the distribution of the virtual axial dipole moment 
are qualitatively well represented in the model. Furthermore the statistics of the times 
between reversals is similar in the model and in the case of the Earth's magnetic field. 
The mean time between reversals depends sensitively on the model parameters. Thus the 
drastic changes in the reversal frequency of the geomagnetic field could be explained by a 
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Figure 11: Top (a): Variability of the virtual axial dipole moment (VADM) during the past 
2 Myr, the SINT-2000 data set. The absolute values of the VADM are given, disregarding 
the five reversals during this period. Bottom left (b) : Power spectrum of the VADM time 
series. Bottom right (c): Distribution of the VADM. 



moderate change of the model parameters with time. 

Our model provides a convincing statistical representation of the geomagnetic field 
reversals process. One should be careful, however, when interpreting the model properties 
in terms of magnetohydrodynamics. Secular variation, which is mainly determined by the 
details of the convective flow dynamics, is certainly not captured correctly. The view that 
the convective columns to a certain degree represent building blocks of the full dynamo 
process seems to be strengthened by our results. A stable polarity can only established 
when the majority of these entities cooperate and produce field of the same polarity. 
Random forcing counteracts this and may sometimes be violent enough to cause a spin to 
flip significantly and leave the team. This may cause its neighbours to follow and ultimately 
lead to a reversal. The magnetic upwellings identified in full 3D dynamo simulations by 
Aubert et al. (2008) could be these events. When these upwellings last long enough or 
produce enough inverse field, they disrupt the normal dynamo process. The statistics of 
the complex interplay of many agents seems to be nicely describable by our domino model 
of Vaks-Larkin type of a set of interacting magnetic spins. 
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